% Copyright (C) 2009 Indrek Mandre <indrek(at)mare.ee>

% usage: [SAMPLES, MINS] = lsexp (A, X, COVM, RANGE, INTERVALS, OPTIONS)
%
% Find the best exponent D corresponding to the equations:
%   A(1)**D =~ X(1,1)    ...    A(1)**D =~ X(1,M)
%   ...
%   A(N)**D =~ X(N,1)    ...    A(N)**D =~ X(N,M)
% This is done by sampling the given range for D using
% the least square minimized residuals.
%
% Arguments:
%   A - N x 1 column matrix of base rows
%   X - N x M matrix of data rows
%   COVM - M x M x N covariance matrixes for the M data elements
%   RANGE - array [low high] specifying the searched range
%   INTERVALS - int, sample by dividing the range by this
%   OPTIONS - optional options structure
%     OPTIONS.info - true/false, print out additional info
%
% Returns:
%   SAMPLES - a matrix with each row containing the sampled
%     dimension and the corresponding residual
%   MINS - 4 column matrix of best exponents. The columns are:
%     dimension - the value of minimums in the range
%     residuals - the minimization value at the given dimension
%     error - error of the dimension using 95% chi-square confidence
%     goodness - 1 if pure chi-square confidence could be used
%       0 if some dubious approximation was used
%
% Required Octave-Forge packages: statistics, optim
%

function [SAMPLES, MINS] = lsexp (A, X, COVM, RANGE, INTERVALS, OPTIONS)
  dof = size(X,1) - size(X, 2) - 1;
  if size(COVM,1) ~= size(COVM, 2) || size(COVM, 3) ~= size(X, 1) || size(A, 1) ~= size(X, 1)
    error('lsexp: Invalid arguments A, X or COVM!');
  end
  if dof <= 0
    error ('lsexp: Not enough rows of data available!');
  end
  chi2tail = chi2inv (0.95, size(X,1) - size(X, 2) - 1);
  SAMPLES = [];
  MINS = [];
  step = (RANGE(2) - RANGE(1)) / INTERVALS;
  for d=RANGE(1):step:RANGE(2)
    S = lsexp_min(d, A, X, COVM);
    SAMPLES = cat (1, SAMPLES, [d S]);
    % if there is a minimum here, find its exact spot
    if size(SAMPLES, 1) >= 3 && SAMPLES(end - 2, 2) > SAMPLES(end - 1, 2) && SAMPLES(end, 2) > SAMPLES(end - 1, 2)
      % use the golden section search to find the minimum x
      [x, y] = fminbnd (@(x) lsexp_min(x, A, X, COVM), SAMPLES(end - 2, 1), SAMPLES(end, 1));
      % find the error using chi2 distribution if possible
      if y < chi2tail
        searchy = y + chi2tail;
        good = true;
      else
        searchy = y + chi2tail;
        good = false;
      end
      % find proper stepsize
      ss = step / 4;
      c = 4;
      while c > 1
        l2 = lsexp_min(x - 2 * ss, A, X, COVM);
        l1 = lsexp_min(x - ss, A, X, COVM);
        r2 = lsexp_min(x + 2 * ss, A, X, COVM);
        r1 = lsexp_min(x + ss, A, X, COVM);
        if l2 > l1 && l1 > y && r2 > r1 && r1 > y
          break
        end
        ss = step / 2;
        c = c - 1;
      end
      [xoff1, tmp, info1] = fzero(@(x) lsexp_min(x, A, X, COVM) - searchy, x + ss, struct ('abstol', 1e-8));
      [xoff2, tmp, info2] = fzero(@(x) lsexp_min(x, A, X, COVM) - searchy, x - ss, struct ('abstol', 1e-8));
      %fprintf (1, '1: at %f: %.8f (%.8f)\n', xoff1, lsexp_min(xoff1, A, X, COVM), searchy);
      %fprintf (1, '2: at %f: %.8f (%.8f)\n', xoff2, lsexp_min(xoff2, A, X, COVM), searchy);
      abserr = max(abs(x - xoff1), abs(x - xoff2));
      MINS = cat(1, MINS, [x y abserr good]);
      if OPTIONS.info
        fprintf (1, '[%7.4f, %7.4f]: %11.8f (y=%.3e), err=%.8f, g=%d\n', SAMPLES(end - 2, 1), SAMPLES(end, 1), x, y, abserr, good);
      end
    end
  end
  SAMPLES = cat (1, SAMPLES, MINS(:, 1:2));
  [tmp, idx] = sort(SAMPLES(:, 1));
  SAMPLES = SAMPLES(idx, :);
end

